

** make bootstrap cis
/*
1. append together bootstrap sample
2. in each sample-week, calcualte bounds based on non-icli and clear cause hospitalizations
3. find 2.5th and 97.5th percentile of bounds (by week) 
4. find average lb and ub, by week, and get weekly bias in bound 
*/

** 1. append together
local B 500
clear
forvalues bs = 1/`B' {
	append using bootstrap_ests/bs_`bs'
}

** 2. calculate intersection bounds 
gen lb_non_m = lower_bound1
gen lb_clear_m = lower_bound1
gen lb_icli_m = lower_bound1 


gen lb_non_i = max( lower_bound1, lower_bound2)
gen lb_clear_i = max(lower_bound1, lower_bound3)
gen lb_icli_i = max(lower_bound1 , lower_bound4)

gen ub_non_m = min(upper_bound1, upper_bound2)
gen ub_clear_m = min(upper_bound1, upper_bound3)
gen ub_icli_m = min(upper_bound1, upper_bound4)

gen ub_non_i = ub_non_m 
gen ub_clear_i = ub_clear_m  
gen ub_icli_i = ub_icli_m 

gen ub_pool_m = min(upper_bound1, upper_bound5)
gen ub_pool_i = ub_pool_m
gen lb_pool_m = lower_bound1
gen lb_pool_i = max(lower_bound1, lower_bound5)

** 3. find percentiles
local B 500
local p025 = ceil(`B'*.025)
local p975 = ceil(`B'*.975)
foreach t in non clear icli pool {
	foreach suff in i m {
		sort week lb_`t'_`suff'
		by week: gen p025_`t'_`suff' = lb_`t'_`suff'[`p025']
		sort week ub_`t'_`suff'
		by week: gen p975_`t'_`suff' = ub_`t'_`suff'[`p975']
	}	
}
sort week lower_bound1
by week: gen p025_pop = lower_bound1[`p025']
sort week upper_bound1
by week: gen p975_pop = upper_bound1[`p975']

sum upper_bound1 if week == mdy(12,11,2020), d
sum p975_pop if week == mdy(12,11,2020)


** find means
foreach t in non clear icli {
	gen flip_`t' = ub_`t'_m == upper_bound1
}
sum flip_* 

sort week 
collapse (mean) lb_non_m-lb_pool_i flip* p025_* p975_*, by(week) 


foreach t in non clear icli {
	foreach b in lb ub {
		foreach suff in i m {
			rename `b'_`t'_`suff' mean_`b'_`t'_`suff'
		}
	}
}
save tables/bootstrap_ci_bias, replace

** set up point estimates and cis   estimates 
use bootstrap_ests/bs_0 , clear
format week %td
format lower_* upper_* %5.4f

gen lb_non_m = lower_bound1
gen lb_clear_m = lower_bound1
gen lb_icli_m = lower_bound1 


gen lb_non_i = max( lower_bound1, lower_bound2)
gen lb_clear_i = max(lower_bound1, lower_bound3)
gen lb_icli_i = max(lower_bound1 , lower_bound4)

gen ub_non_m = min(upper_bound1, upper_bound2)
gen ub_clear_m = min(upper_bound1, upper_bound3)
gen ub_icli_m = min(upper_bound1, upper_bound4)

gen ub_non_i = ub_non_m 
gen ub_clear_i = ub_clear_m  
gen ub_icli_i = ub_icli_m 


gen ub_pool_m = min(upper_bound1, upper_bound5)
gen ub_pool_i = ub_pool_m
gen lb_pool_m = lower_bound1
gen lb_pool_i = max(lower_bound1, lower_bound5)


rename lower_bound1 lb_pop
rename upper_bound1 ub_pop

merge 1:1 week using tables/bootstrap_ci_bias, assert(3) 
drop _merge

** bias ub m 
foreach t in non clear icli {
	gen bias_ub_`t' = mean_ub_`t'_m - ub_`t'_m 
	gen bias_lb_`t' = mean_ub_`t'_i - lb_`t'_i
}


format week %td 

save tables/estimates_bs_ci, replace 

